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method for delay optimal control 
problems 
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Abstract 


In this paper, we present an efficient method to solve linear time-delay 
optimal control problems with a quadratic cost function. In this regard, 
first, by employing the Pontryagin maximum principle to time-delay sys- 
tems, the original problem is converted into a sequence of two-point bound- 
ary value problems (TPBVPs) that have both advance and delay terms. 
Then, using the continuous Runge-Kutta (CRK) method, the resulting 
sequences are recursively solved by the shooting method to obtain an opti- 
mal control law. This obtained optimal control consists of a linear feedback 
term, which is obtained by solving a Riccati matrix differential equation, 
and a forward term, which is an infinite sum of adjoint vectors, that can 
be obtained by solving sequences of delay TPBVPs by the shooting CRK 
method. Finally, numerical results and their comparison with other avail- 
able results illustrate the high accuracy and efficiency of our proposed 
method. 
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1 Introduction 


In recent years, optimization and control of systems with time delay have 
been considered in much research because the time delay in many processes 
cannot be ignored. To more accurately express the behavior of a natural 
phenomenon, we need a more complex system. Some of the applications of 
these issues are in the chemical, electronic, medicine, engineering, biological, 
economy, and so on [22, 19, 12, 7, 8, 41]. 

In general, two methods are provided to solve optimal control problems 
(OCPs). The first approach involves the use of necessary and (or) sufficient 
conditions of optimality by applying the Pontryagin minimum (maximum) 
principle or optimality principle. The minimum principle was presented in 
1956 by the Russian mathematician Lev Pontryagin and his students, and 
its primary application was to maximize the terminal velocity of a rocket. 
This result was obtained using the classical ideas of variational calculus. 
The equations obtained from these conditions can be solved numerically. 
This approach yields indirect methods, which are known as analytical-based 
methods; see [39, 43, 17, 13]. 

In another approach, an OCP is considered an optimization problem. in- 
stead of using the optimality conditions, the dynamic constraints are trans- 
formed into an algebraic equations system by discretizing the time interval 
and parameterizing the variables of the problem. Therefore, the OCP be- 
comes a nonlinear programming problem of dimension finite. The result- 
ing nonlinear programming problem can then be solved using optimization 
techniques. This approach yields direct methods. We refer the reader to 
[11, 2, 18, 26, 8]. Since direct methods do not need to calculate the opti- 
mality conditions, they can be used for a wide range of OCPs. However, 
the lack of guarantee for the optimal solution and the high amount of mem- 
ory resources and time for producing a close approximation is among the 
disadvantages of these methods. 

In the case of time-delay OCPs, in 1963, Oguztéreli [35] was one of the 
pioneers in the analytical-based approach (also, see [36]). For the first time, 
Kharatishvili [24] generalized the Pontryagin maximum principle for OCPs 
with a constant delay in the state variable. Then in [25], he gave similar 
results on OCPs with delay in the control variable. After that, in 1968, a 
maximum principle for OCPs with multiple constant delays in state and con- 
trol was proved by Halanay [16]. In 1972, Ray and Soliman [42] also obtained 
similar results. Guinn [14] transformed the delayed OCP with constant delay 
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in the state variable into a higher-dimensional undelayed OCP. Banks [3] de- 
rives a maximum principle for control systems with a time-dependent delay 
in the state variable. 


The system resulted from the necessary conditions that Kharatishvili 
provided, which was a two-point boundary value problem involving both 
advance and delay terms. This type of problem does not have an exact 
solution, except in exceptional cases. Therefore, there are many attempts 
available in the literature to approximately solve this problem; for example, 
see [29, 44, 30, 31, 32, 20, 21, 6]. 


The following articles can be mentioned as the latest studies. For OCPs 
with time-invariant delayed systems, Mirhosseini-Alizamini, the second au- 
thor, and Heydari [32] applied the variational iteration method and then 
obtained a suboptimal solution for the two-point boundary value problem 
(TPBVP). Moreover, Mirhosseini-Alizamini and the second author [31] in- 
vestigated infinite horizon OCPs with time-variant delayed systems. Also, 
using a Hermite interpolation polynomial for delay terms and employing a 
second-order finite difference formula for the first-order derivatives, Jajarmi 
and Hajipour [21] converted the TPBVP obtained from the time-delay OCP 
into a system of linear algebraic equations and then solved it. Recently, using 
an algorithm based on the forward and backward difference approximation, 
Bouajaji et al. [6] solved the system obtained from the application of the 
Pontryagin maximum principle to a delayed OCP. 


In this work, we investigate a family of time-delay OCPs with a quadratic 
cost functional that should be minimized subject to a linear time-delay system 
with constant delay in the state variable. Using the Pontryagin minimum 
principle for delayed systems from [24] and then applying continuous Runge— 
Kutta (CRK) methods, we convert a time-delay OCP into a sequence of 
linear TPBVPs and thereafter solve it recursively by the shooting method to 
obtain the optimal control law. 


The rest of the paper is organized as follows: The CRK methods are pre- 
sented in Section 2. After that, in Section 3, we introduce the Shooting CRK 
(SCRK) method and apply it to a delayed TPBVP. Then, in the continua- 
tion of this section, we present a basic algorithm for the proposed method. 
In the next section, we will use a generalization of this algorithm to solve a 
time-delay OCP. Section 4 describes the Pontryagin maximum principle for 
our delayed OCP and designs an algorithm based on the previous algorithm 
defined in Section 3 for solving the final system. In Section 5, we give sev- 
eral numerical examples to demonstrate the effectiveness and accuracy of the 
proposed technique. Finally, with the conclusion in Section 6, we end the 
article. 
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2 CRK methods 


In this section, we describe the CRK methods. Consider f(t,a(t)) € 
C°([to, tr] x R¢,R*). The CRK methods were originally designed to treat 
the initial value problem for the following ordinary differential equation: 


&(t) = f(t,e(t)), to<t<ty, 
_ (1) 
x(to) = 2. 

Some of the implicit Runge-Kutta methods are equivalent to collocation 
methods; see [46]. Thus, they sequentially provide a continuous extension 
of the approximate solution without any additional evaluation of f. Indeed, 
the next question is whether there is such a continuous extension for each 
Runge-Kutta process that is given sequentially by the method itself? 

Nersett and Wanner partially answered this question by proving that a 
large number of Runge-Kutta methods are the same as the somewhat per- 
tubed collocation method that is somewhat perturbed. After that, Zennaro 
[47] presented a continuous extension of the solution provided by a Runge— 
Kutta method, which includes the collocation solution if it is equivalent to 
collocation and behaves similarly in other cases. 

Let A = {to,...,tn,...,tw = ty} be an arbitrary mesh. Then for the nu- 
merical solution of the ordinary differential equation (1), an s-stage discrete 
Runge-Kutta method has the form 


s 
Tr41 = In + hn41 SS bik, (2) 
i=l 
8 
ky = ftw haa >, Gh); i= | rere P (3) 
g=l 
where 7 oi a;;, tt, = tn + CiAn4i1,t = 1s; +58, and An4t — tn41 == bins 


In addition, the Runge-Kutta method (2) and (3) is denoted by (A, 6). Let 
the solution have advanced to the point t = t,. Zennaro [47] showed that 
for this s-stage Runge-Kutta method of order p, there is a CRK method of 
degree d, if there exist s polynomials 6;(@), 7 = 1,...,s, of degree less than 
or equal to d, independent of f. This method reads as follows: 


n(tn + Ohngi) =2n +hnsi >> bi(A)ki, O<O<1, (4) 
t=1 
ke = f (thy an + >— aigky), t= 12g 8} (5) 
j=l 


where 
n(tn) =n; n(tn am An+1) = In+1; 
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and 2, is an approximate solution obtained by applying the Runge-Kutta 
method for x(t,). This method, which is usually expressed as (A, b(0)), can 
also be related to the following CRK tableau: 


c| A 
oT (8) 


In fact, {c;, ai; }’s are the same as the coefficients of the discrete Runge-Kutta 
method. Now, we recall the consistency of the discrete Runge-Kutta method 
from [5]. 


Definition 1. [5, Definition 5.1.3] Consider p > 1 the largest integer having 
the following property: For every mesh point and C?-continuous right-hand- 
side function f(¢,2) in (1), the local solution z,+1() to the local problem 


tier - Ftivzntalt) bi StS as re 


n? 


satisfies 


1 
eee (he) —eagal| = On) 


uniformly with respect to 2* belonging to a bounded subset of R@ and re- 
spect ton = 0,..., N—1, Then we say that the discrete Runge-Kutta method 
(A, b) is consistent with order p. 


Similarly, with the above notations, we say that the continuous extension 
(4) is consistent with uniform order q if g > 1 is the largest integer having 
the following property: 


41 
go ea) — n(t)|| = O(hA), 


for every mesh point and C%-continuous right-hand-side function f(t, x) in 


(1). 


According to Definition 1, the convergence results in discrete and CRK 
methods for ordinary differential equations have been proved in the following 
theorem; see [5]. 


Theorem 1. [5, Theorem 5.1.4] Suppose that the Runge-Kutta method (2) 
and (3) is consistent with order p and that f(t,2) defined in (1) is a right- 
hand-side C?-continuous function. Then, on any bounded interval [to, t;], 
the method has discrete global order (or, equivalently, is convergent of order) 


p. In other words, 


= = Pp 
max, len) Lp|| = O(h?), 
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in which h= max hy. 
1<n<N 
Moreover, let the continuous extension (4) have the uniform order g. Then 


the CRK method (4) and (5) has the uniform global order (or, equivalently, 
uniformly convergent of order) g’ = min(q+1,p), which means that 


max ||x(t) — n(t)|| = O(h7). 


toStcty 


Then, Baker and Paul [1] generalized this idea for a CRK method to delay 
differential equations with a general delay differential equation of the form 


£(t) = ft, 2), 2 —7()), t> to, (7) 
z(t) = d(t), to — T(to) < t < to, 

in which f : Rx R" x R” > R” and r(t) > 0. Moreover, ¢ € O° [to —T(to), to] 

denotes the initial information of the state variable x. For delay differential 
equations, Baker and Paul [1] modified (4) and (5) as follows : 


n(tn + Ohnti) =2n+hnti > di(O)ki, OS O<1, (8) 
i=l 

X; = In + hms > aigh;, $= Lys day 8! (10) 
j=l 


When the delay is constant and hy41 <7, then 7(t?, — rT) is known for any 7 
(0 < ¢ <1). In this case, 7(t!, — 7) is available from the past, so this method 
is an explicit CRK method. The pair formed by (A, b) and (A, b(@)) is called 
the underlying CRK method. 


Theorem 2. [5, Theorem 6.3.1] Assuming the delay differential equation (7), 
suppose that f(t,z,y) € [to,t7] x R” x R” is a C?-continuous function. Then 
the delay 7(t) € [to, t7] x R” is a C?-continuous function and ¢(t) is the initial 
C?-continuous function. In addition, let A = {to,ti,...,tn,...,tw =ty} be 
the mesh containing all points of discontinuity with the order less than or 
equal to p being in [to,ty]. Also, assume that the underlying CRK method 
has the uniform and discrete orders g and p, respectively. Then for the delay 
differential equation, the CRK method (8), (9), and (10) has uniform global 
and discrete global orders q’ = min (p,q +1). In other words, 


lea l= Oh), 


and 
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= =x qd’ 
13x, l(b) — nll = OCR"), 


where h= max hy. 
1<n<N 


3 Outline of SCRK method for a delay TPBVP 


In the present section, we first state details of the proposed method on a 
TPBVP with only a time-delay term. Therefore, consider the following basic 
form of a first-order TPBVP with a time delay: 


z(t) = filt, c(t), y(t), c(t —7), t—ir)), to <t< ty, 

y(t) = fo(t, x(t y(t pe Se eyes 7))s to St Sty, (11) 
x(t) = d(t), to-—TS<t< bo, 

y(tr) = 6 


For solving this problem, we need to use the solutions to a sequence of ini- 
tial value problems that are made by substituting the initial guess y(to) = z 
instead of the terminal condition y(t) = 8 in (11). 


To approximate a solution to the boundary value problem (11), we involve 
a parameter z, by choosing the parameters z = zz such that 


jim y (ty, zn) = y(ts) = B, 


where y(t) is the solution to the boundary value problem (11) and y(t, zx) 
denotes the solutions to the constructed initial value problem with initial 
conditions x(t) = o(t), to -T <t < to and y(to) = zx. 

This technique is called the Shooting method. For starting, we choose 
a parameter z, such that it determines the initial evaluation at which the 
object is fired from the point (to, d(to)) and along the curve indicated by the 
solution to the problem 


a(t) = filt, c(t), yt), 2(t— Jult —7)), to St <ty, 

y(t) = fo(t, x(t Y t ,clt= ),ylé—7)), to St< ty, (12) 
a(t) = o(t), ty -T<t< to, 

y(to) = a 


If y(tp, 21) is not sufficiently close to 8, then we correct the approximation 
by choosing elevations z2, z3, and so on, until y(t;, z,) is sufficiently close to 
B. 


For determining the parameters z,, we must solve this problem: 


y(t, 2) — B= 0. (13) 
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To solve this nonlinear equation, we use the secant method. For this method, 
we need to choose initial approximations z, and z2 and then generate the 
remaining terms of the sequence by the following formula: 


y(te, Zk) aa B 


Zk —Zk-1), K=3,4,.... 14 
ieee a4) 


Zkt+1 = 2k 


To obtain y(t, 21) in (12), we use the CRK method (8), (9), and (10) for a sys- 
tem of delay differential equations. For a given mesh A = {to,...,tn,...,tnw = 
tr}, leth= ua. In each underlying mesh interval [t,,tn41], CRK formulas 
for (12) are as follows: 


ki = Fits Magee ee et, Ti) Hilby oa T2)), i= 1, oe 5, 
kai = fo(tt,, Xi, Yi, Ne (t, a, 71) My (th, _ T2)), t= 1, oe 5, 
Xi = In +h, aigh,j, ol ee 


. ; (15) 
Yi = yn thins Gigko,j, join Sena 2 
Ny (tn + Oh) = Yn +h dy_, bi(O)ke 3, 0<é<1, 


Note that at the endpoint of the interval, the stop condition must be 
checked. In the following algorithm, we describe an SCRK method for a 
time-delay TPBVP. 


Example 1. Consider the following second-order delay boundary value prob- 
lem: 


uw 


x (t) =—#sina(t)-— (t+ 1l)e(t-1) +t, 0<t<2, 
a(t)=t—$, t<0, (16) 
x(2) = —t. 


With the new condition x (0) = z, instead of solving (16), we need to solve 
a sequence of initial value problems of the form 


x (f) =—#sine(t)-— (t+ 1l)et-1 +t, O<t <2, 
a(t) =t— 5, <0; (17) 
x (0) =z 

Now, we try to make the value of y(2, z) as close to 8 = —3 as possible by 


adjusting the value of z. Before that, by assuming y(t) = a (t), we turn the 
delay second-order system (17) into a delay first-order system as follows: 


x (t) = y(t), oreo: 
y (t) = —% sin a(é) — (t+ 1)x(t- 1) +t, 
y(0) =z 
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Algorithm 1 SCRK method for time-delay TPBVP 


Step 1. 


Step 2. 


Step 3. 


Step 4. 


Step 5. 


Step 6. 


Step 7. 


Step 8. 


Set N (the number of subintervals), h = ‘t0 K =1,M (the maximum number 
of iterations), and s (the number of stages of the CRK method), and choose 21, z2, 
and tolerance error bound e. 

While L < M, do 


e Set ro =a and y(to) = 21, 
For k = 1,2,---, 
e solve (12), using the CRK method (15). 
e Set ro =a and y(to) = 21, 
Check the stop condition, 
e If lynx — B| < «, then the procedure is complete, and jump to Step 7, 
e else, go to the next step. 
If k = 1, then set y(to) = z2 and back to Step 3, 
e else, go to the next step, 


Calculate the next approximation for z,4, from (14), set y(to) = z,41, and back 
to Step 3. 


e end for 
Stop the algorithm and output (tn, ¢n, yn). 
e end while 
Output (maximum number of iterations exceeded). 


e Stop 


We solve this problem by applying Algorithm 1. For this purpose, we use the 
explicit Runge-Kutta of discrete order p = 4 with the following coefficients: 


0 | 0 
i]2 ° 
2 0 5 O 
5 [0-0 1.0 
| re ee es 
6 3 3 6 
Moreover, we set 
1 2 1 
b1(0) = 07 + =6 b3(0) = =0 
i( ) 5) + 3 ’ 3( ) 3 ’ 
1 1 
bo(0) = 0 ba(0) = ~0? — <0. 
2 ( ) 3 ’ 4( ) 3 
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Table 1 indicates a comparison between the approximate result of our 
SCRK method and the results obtained in [34]. 


Table 1: Approximation values of x(t) in Example 1 


By (1) Ly (1.5) En (2) 
n Ref. [34] Proposed method Ref. [34] Proposed method Ref. [34] Proposed method 
4 -1.854384 -1.983957 -1.719174 -1.884111 -0.499976 -0.499999 
6 -2.018854 -2.032385 -1.896332 -1.922809 -0.499999 -0.500000 
8 -2.066385 -2.052802 -1.946231 -1.939199 -0.499999 -0.500000 
10 -2.078723 -2.063029 -1.959110 -1.947469 -0.499999 -0.500000 
12 -2.081821 -2.068830 -1.962343 -1.952141 -0.500000 -0.500000 


4 Design of SCRK method for an OCP with time delay 
in state variable 


In this section, we first use the Pontryagin maximum principle to solve our 
delayed OCP. Then, for solving the final system, we describe an algorithm 
based on Algorithm 1. Through this section, by PC'([to,ts],R”) we de- 
note the class of continuous functions from [to,t;] into R” whose first-order 
derivatives are piecewise continuous, and similarly, PC((to,ts],R”) denotes 
the class of piecewise continuous functions from [to, ty] into R”. 


Consider the linear system with delay in the state variable 


Ax(t)+ Aia(t—7)+ Bult), to<t<ty, 


(19) 
to -TStS bo, 


———N 
8 2 
—_~— 
cH oe 
eset: Seal: 
ol 
a 
e 
=o 


where u(t) in PC([to, tr], R”) and x(t) in PC1([to —7, ty], R”) are the control 
and state variables, respectively. In fact, the parameter 7 > 0 is nonnegative 
and indicates the time delay. Furthermore, the initial state function ¢(¢) is 
continuous in C([to — 7, to], R”), and finally, the matrices A, B, and A, are 
real constants with appropriate dimensions. For t € [to,t,/], our aim is to 
obtain, u*(t), the optimal control law minimizing the quadratic cost function 


1 


Jaz f WP @Rule) +2" OQeOa+ 5e™ENQsA(t)), CD) 


in which R € R™*” is a positive definite matrix and Q and Qy € R"*” are 
positive semi-definite matrices. 


For time-delay OCPs, it follows from [24] that the pontryagin maximum 
principle provides the necessary conditions of optimality for the problem (19) 
and (20) as follows: 
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a(t) = Aa(t) + Ayx(t —T) — BR-1BT X(d), 16 SEAT 
ox ees a ee 

—Qzx(t) — AT Ad), tp—T<t<ty, 
a(t) = 9(t), to -T <t< to, 


Ate) = Qru(ty). 
(21) 
The Hamiltonian function from which the above conditions are derived is 


H(a,u,d,t) = 7 (t)[Ax(t) + Bu(t) + Aia(t — 7) + 527 (Qz(t)+ 
sit (t) Ru(t)), (22) 


where A(t) € PC1([to, ty], R”) is called co-state vector. Moreover, 
u*(t) = —R7'BT X(t), (23) 


for to < t < ty, is the optimal control law. We recall that the system (21) 
is a TPBVP with both time-advance and time-delay terms. Unfortunately, 
in general, this problem does not have any analytical solution. Therefore, 
providing an efficient method for solving this difficult problem numerically is 
very important. 


At first, we produce a sequence of TPBVP as 


é™)(t) = —SA) (¢t) + Ax (t) + Arx™(t—7), to St < ty, 
4, = ae — Qe (t)— ATAR VET), toStStp—7, 
—AT \®(t) — Qz™ (£), tp -—T<tK<ty, 
x) (t) = 9(t), to-T<t< to, 
AM (Ep) = Qa (Ey), 
2OW=0, ADM =0, tr <t<t,, 
(24) 
where S = BR~!B? and k = 1,2,.... Therefore, 


u*) (t) = —R-1 BT) (4) (25) 


is the sequence of controls. Now, we are ready to obtain a closed-loop optimal 
control. We can define the co-state vector by 


AM) (t) = g) (4) + P(t)z™ (2), (26) 
in which g(t) € R” is the kth adjoint vector and P(t) € R”*” is an 


unknown function matrix with positive-definite property [45, 44]. 
Consider the following extended sequence of the TPBVP (24): 
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£)(t) = [A — SP(t)|2™ (t) — Sg (t) + Aix (t — 7), to <t<ty, 
—P(t)Aia™ (t — 7) — [A— SP(t)]7 9 (t) 
g(t) = ¢ —APP(t+ r)xX@-Y(t+7) — ATg@-U(t+7), to <t<ty—7, 
—P(t)A;z™ (t — 7) — [A— SP(t)]7 9 (2), tp -—T<t<ty, 
c(*)(t) = g(t), fo -T<tX< to, 
g(t) = 0, 
c(t) =0, g(t) =0, to <t<ty. 
(27) 


We note that by substituting (26) into the first equation of (24), the kth 
optimal closed-loop system is constructed, which is the first equation of the 
system (27). Similarly, substituting (26) in the second equation of (24) and 
comparing the result with the derivative of (26), we obtain the second equa- 
tion of the system (27). Also, 


—P(t) = P(t)A+ A’ P(t) — P(t)BR-'B? P(t) +Q, 
P(tr) = Qs, (28) 


is a Riccati matrix differential equation. 
Moreover, from (25) and (26), the sequence of controls is converted to 


u®) (¢) = —R1 BT (Pc (t) + 9(t)), k=1,2,.... (29) 


The system (27) is similar to (11), except that (27) has advance terms in 
addition to the delay terms. Now, we want to use Algorithm 1 to solve this 
advance-delay TPBVP. By using the SCRK method, we have the following 
CRK iteration formula of (27) in the mesh interval [t,,, tr+1]: 


1 (tn + Oh) =a +>” :(O)[WEE (ah) + b>” aighs,s) 
i=1 j=l 


— Sg) +h>° aizko,s) + Ain (t, —7)], to St < ty, 
j=l 
(30) 
(E) 1 ays. (6)|—-WT 42 )(gQ 4+hyve oy ae 
gn +hYT-, bi(4)L (tr, )(gn > + jai Vs 2,3) 
—P(t,) Aine” (th, — 7) — AT P(t, + 7)2®-Y i, +7) 


(tr +Oh)=4 —ATg@ VL +7), to St<ty—7, 
k s 4 k s 
gh? + RIF b:(8)[— WT (ti) (gh? + hDOS_, aizho,s) 
—P(ti,) Ain (ti, — r)], tp —T <t<ty, 


(31) 


where t?, = tn +cih, U(t?,) = A— SP(tt,), and 0 < @< 1. Also, 
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w(t) = $(t), to—T StS to, (32) 
g(t 7) = 0; 
are the known initial and final conditions. 

As already mentioned, for the constant delay, if 0 <c; < 1 and h <1, 
then n(t?, —7) is known for any i. Hence, 7,(t?,—7) in (30) and (31) is known, 
and there is no so-called overlapping. On the other hand, «‘*~!) (t?, +7) and 

(kV) (¢?, 4+ 7) are obtained from the previous iteration by the assumptions 
2 (¢t) = 0 and g(t) =0. 


Theorem 3. Consider TPBVP (27). 


i) Assume that the right-hand-side functions corresponding to #")(t) and 
g(t) and the initial function ¢(t) are C?-continuous in their domains 
(p is the discrete order of the underlying CRK method). Then the 
sequences {n*) (t)} and {ni*) (t)} obtained from CRK formulas (30) 
and (31) with initial and boundary conditions (32), converge uniformly 
to the solution of TPBVP (27). 


ii) Under the assumptions of part (i), the sequences {u*)(t)} and {J}, 
which are defined as follows 


u(t) = —R*B™ [P(t)n (t) +n ()], (33) 
le sn tg) pnt (ty) © ; Ta) Qn® 0) 
+ (u(t)? Ru (t)] dt, (34) 


converge to optimal control u*(t) and the optimal value of objective 
function, J*, respectively. 


Proof. i) Consider the vector function F' as follows: 
F(t, 0,9, u,v, w, 2) = (&(t), g(t)", to St < ty, 


and u,v, w, z denote delay and advance terms corresponding to the vari- 
ables x(t) and g(t). Also, «(t) and g(t) are the functions defined in (27). 
Because F and ¢ are C?-continuous functions and 7 is a constant delay, 
according to Theorem 2, the sequences {nh*) (t)} and {ni*) (t)} from the 
CRK method are uniformly convergence to the exact solutions of (27). 


Suppose that {n\* (t)} and {ns*) (t)} are solution sequences produced 
by the CRK method, which are convergence to #,(t) and #,(t) under 
the assumptions of part (i). We take the limit from the (33) as k + ov, 


Me 
mE: 
—/ 


a(t) = Jim w(t) = —R“BTP(A)( dim nf (1) + Jim nf? (2) 


k—-o0o 
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= —R-'BT[P(t)fie(t) + fg(O). 


Since #,(t) and 7,(t) are the exact solutions of necessary conditions 
(27), so ti(t) is the optimal control u*(t). 
Similarly, we take the limit from the (34) as follows: 

J := lim J) 


k- 00 


= Jim (GOP ep) Ort) 


k-0o 


+5 Jim Cf “(nT Onl) + (u® (O)? Ru® (tat 


k-o0o 
Petes 
=5( jim (nf (t5))7 Qp( lim nf" (t7)) 


pe i "(tim (n®(#))7)Q@(tim (8) 


k-> oo k—-00 
+( 


sim (ul ()T)RC jim. ul (0) 
lor - Lf ramos i” (t) Rai(t)|dt 
= SAE NQ relly) +5 [EQ + AH RAC 


so, J is the optimal value of the performance index J. 


According to Theorem 3, it can be concluded that for enough iterations 
of the CRK method, for example, N iterations, where N depends on a given 
error criterion, we can obtain a suboptimal control as follows: 


u(t) = —R*BT[P(t)n? (t) +n (0)]. (35) 


In this case, the continuous suboptimal state function is as 7,(t) & ns? (t). 
To calculate a more accurate state function, the suboptimal control function 
resulting from equation (35), can be placed in (19), and we then solve the 
obtained initial value problem. Finally, by placing this pair of suboptimal 
control and state in the objective function, we have 


d sie a 
I) = SUE) + 5 f CMO) AN 
0 
+ (u) (t))P Ru (t)) dt. (36) 
For given ¢ > 0, if the stop condition, 


JN) — JON-1) 


si —| <5 
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is satisfied, then the suboptimal control (35) will have the desired accuracy. 
Now, to implement the above method, we provide the following simple algo- 
rithm. 


Algorithm 2 SCRK method for time-delay OCPs 


Step 1. Solve P(t) from (28). 

Step 2. Put k = 1, «© =0, and g) = 0. Then obtain a continuous approximation for 
x(*)(t) and g(*)(t) from the kth TPBVP (30), (31), and (32) with the shooting 
method (Algorithm 1). 

Step 3. Let N =k and obtain u)(t) from (35). 

Step 4. Obtain J(%) from (36). 


JON) g(N-1) 
Step 5. If —a 


e else, let k:=k-+1, and back to Step 2. 


<e, then the procedure is complete, and go to the next step; 


Step 6. Stop the algorithm and consider the output u(%)(t) as the desired closed-loop 
suboptimal control law. 


5 Numerical examples 


Now, we are ready to present several examples for showing the efficiency of 
the SCRK method. 


Example 2. Consider the delay system 


t=a(t)+u(t)+a2(t-1), t>0, (37) 
a(t) =1, =1<4< 0, 
to minimize this quadratic cost functional 
3.5 a ae 
J==a2°(2)+— ] u*(t)dt. (38) 
2 as, 
It follows from [4] that the exact solution for u(t) is 
7 6(e2-*+(1—-tet), 0<t<1, 
a= im teres, oe 


and that J* = 3.1017, where 6 = —0.3932. According to (37) and (38), we 
have Q=0, R=1, Q7 =3, A=1, B=1, and A; =1. Hence, (28) can be 


rewritten as 
p(t) + 2p(t) — p(t) = 0, 
p(2) = 3, 


which has the unique solution 


(40) 
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6e4-2¢ 


= 3230-3" (41) 


p(t) 


For the first time, Banks and Burns [4] proposed a numerical method to 
solve this problem based on averaging approximations. Then Pananisamy 
and Rao [38] solved it by using the Walsh functions. After that, Mirhosseini- 
Alizamini, the second author, and Heydari [32] used the variational iteration 
method. Furthermore, Jajarmi and Hajipour [20] employed a finite difference 
method for solving this problem. We apply our proposed method according 
to Algorithm 2 to this example. Comparison results of the optimal values 
of J obtained by our proposed technique and other mentioned methods are 
listed in Table 2. The curves depicted from the obtained approximations for 
the state and control variables of problems (37) and (38) are shown in Figure 
6. 


Table 2: Value of cost functional for various methods in Example 2 


Method J 
Banks and Burns [4] 3.0833 
Pananismay and Rao [38] 3.0879 
Mirhosseini-Alizamini, Effati, and Heydari [32] 3.1091 
Jajarmi and Hajipour [20] 3.101717 
Proposed SCRK method 3.101667 
Optimal cost J* 3.1017 


0 0,5 1 15 2 


t ——_" Approximation Exact 


(a) State variable x(t) b) Control variable u(t) 


Figure 1: Simulated curves of (a) state variable and (b) approximation and exact values 
of control variable for Example 2 


Now, we give another example. 


Example 3. Consider the time-delay system 
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=u(t)—a(t—-1), 0O<t<l1 
#=u(t)—2(t-1), OS#<1, a 
a(t) =1, S15 630, 
to minimize this quadratic cost functional 
Pail 5 age 
Jie [~a*(t) + ~u*(t)]dt. (43) 
9 2 2 


Now, our aim is to obtain the optimal control, u(t), subject to (42) that 
minimizes (43). Moreover, the Riccati equation for this example is 


p(1) = 0, 


and has the unique solution 
p(t) = —tanh(t — 1) (45) 


The exact solutions for u(t) and x(t) are, respectively, obtained as follows: 


u*(t)=14+ (sinh(t — 1) — cosh(#)), (46) 


il. 
cosh(1) 


x(t) = any ost — 1) —sinh(f)). (47) 

Moreover, it follows from [33] that the optimal value of cost functional 
is J* = 0.1480542786. It can be shown that the approximate value of the 
cost functional calculated by the proposed SCRK method is equal to J = 
0.1480542988. It is clear that the approximate value of J is very close to the 
optimal value. Also, we depict the simulation curves of the trajectory of x(t), 
control variable u(t), and their exact values in Figure 2. 


For the first time, Eller, Aggarwal, and Banks [10] presented the next 
example and then studied by other authors in [23, 37, 9, 40]. 


Example 4. Consider the linear time-varying delay system 


&=a(t)+u(t)+a(t—1), 0<t<2, (48) 
a(t) =1, Ze ai 
to minimize this quadratic functional 
2 
J= | [x?(t) + u?(t)]dt. (49) 
0 


Therefore, the Riccati equation for this example is 
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-0.5 0 0.5 1 


—— Approximation — — Exact 


—— Approximation — — Exact 


(a) State variable x(t) (b) Control variable u(t) 


Approximation and exact values of state and control variables for Example 


Figure 2: 
3 

I(t) + 2p(t) — Zp?(t) +2 =0 

B(t) + 2p(t) — 4p?() +2. =0, i 

p(2) =0, 
and the unique solution for this Riccati equation is 

-1,v2 
p(t) = 2— 2V2tanh(/2¢ + tanh ore ks 2/2). (51) 


In Table 3, we compare the results of the suggested method with the reported 
results in [10, 23, 37, 9, 40, 21]. Figure 3 shows the approximate values of 


the state and control variables of the problem (48) and (49). 


Table 3: Values of cost functional for various methods in Example 4 


Method J 
Eller, Aggarwal, and Banks [10] 6.45 
Dadebo and luus [9] 6.26775 
Oh and Luus [37] 6.23711 
Jamshidi and malek-Zavarei [23] 6.5 
Santos and Sanchez-Diaz [40] 6.97 
Jajarmi and Hajipour [21] 6.219615 
Proposed SCRK method 6.200623 


Example 5. In this example, we want to minimize the cost functional 


(52) 
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(a) State variable x(t) (b) Control variable u(t) 


Figure 3: Simulated curves of (a) state and (b) control variables for Example 4 


with the following two-dimensional delay system: 


&1(t) = x2(t), O<t<2, 
&o(t) = —a1(t) — ao(t-—1)+ u(t), 0<t <2, (53) 
21(0)=10, —9(0) =0, -1<t<0. 


Now, our aim is to obtain the optimal control u*(t) subject to (53) that 
minimizes (52). It follows from [4] that this problem has the exact solution 
dsin(2—t)+ $(1—t)sin(t—1), O<t<1 
ue(t) = [8-H + $A sine), OStEL, 

dsin(2 — t), 1<t<2, 


in which the optimal cost is J* = 3.3991 and 6 = 2.5599. In this two- 
dimensional example, we have A = 1 op A, = F : | B= A Q= 


0-1 
00 100 
9), y= 99), and = 


Thus, instead of the Riccati equation, we have a system consisting of four 
equations and four variables. After applying the proposed method to this 
example, we obtained the minimum value of J = 3.3993. In Table 4, the 
comparison of the result obtained with our proposed method and the result 
based on the techniques presented in [4, 28, 27, 15, 32] is shown. Also, Figures 
3 and 5 show the corresponding state trajectories of x(t), v2(t) and control 
variable u(t), respectively. 
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Table 4: Cost functional values of various methods for Example 5 
Method J 
Banks and Burns [4] 3.2587 
Lee [28] 3.4827 
Khellat [27] 3.43254 
Haddadi, Ordokhani, and Razzaghi [15] 3.21663 
Mirhosseini-Alizamini, Effati, and Heydari [32] 3.3991 
Proposed SCRK method 3.3993 
107 0 
-14 
3 
-24 
6 
x -34 
1 
4 ? -44 
24 5 
-67 
0 + + + 
0.5 1 1.5 2 
t “7 
(a) State variable x1 (t) (b) State variable x2(t) 


Figure 4: Simulated curves of state variables for Example 5 


0 0.5 i 15 2 


—— Approximation — — Exact 


Figure 5: Control variable u(t) for Example 5 


6 Conclusion 


We employed The CRK method to solve a class of time-delay OCPs with 
delay in the state variable and with quadratic cost functional in this paper. At 
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first, by employing the Pontryagin maximum principle for time-delay systems, 
the delay OCP was converted to a sequence of TPBVPs that have both 
delays and advance terms. After that, by applying the CRK method together 
with the shooting method, we constructed two sequences in which the delay 
and advance terms are known. Then we showed that by establishing the 
continuity condition, these sequences converge to the exact solution of the 
problem. The numerical results were presented to illustrate the high accuracy 
and efficiency of our proposed approach. Further research can be done on 
the extension of the SCRK method for solving time-delay OCPs with time- 
dependent delays in the control and state variables. 
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